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^ ■ Abstract. We carry out the Ewald summation for the Rotnc-Prager-Yamakawa 

mobility tensor, the Oseen mobility tensor and further variations of both, relevant for 
the hydrodynamic interactions in colloidal suspensions, where all interacting particles 
are within a single plane, i.e., adsorbed at a fluid interface or other quasi two— 
'■^J | dimensional systems. We use the Poisson summation formula for systems periodic 

■ in two dimensions and finite in the third dimension in order to obtain simple formulae 

for applications, such as molecular dynamics or Brownian dynamics simulations. We 
show, that for such systems, as soon as noise is taken into account, a commonly 
used approximate three-dimensional Ewald summation leads to a spurious system size 
\ dependence, which may considerably affect the interpretation of simulation results and 

' will be cured within our approach. Additionally, the resulting formulae are found to be 

computationally much less expensive than the approximate three-dimensional Ewald 
summation. 
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1. Introduction 

The influence of hydrodynamic interactions (HI) on the dynamics of colloidal systems 
(i.e. colloidal suspensions or colloids trapped at an interface) is subject to ongoing 
research [TJ |2} EJ HJ |5j El [7], both from theoretical and experimental points of view. 
In many circumstances, these many body systems can be investigated only with the 
help of simulations. For colloids floating in a bulk solvent at low Reynolds number, a 
reasonable treatment of HI can be achieved within Stokesian dynamics [8] . In particular 
the Oseen or Rotne-Prager-Yamakawa far-field approximation [9| [TO] . treating the HI as 
pairwise additive interactions, allows an implementation of Hydrodynamic Interactions 
within Brownian dynamics computer simulations [H [TTJ. Since the hydrodynamic 
interactions in the bulk decay oc 1/r, where r is the distance between particles (a 
similar component is also present in the vicinity of or at fluid-fluid-interfaces [32} [13] ) , 
they are considered to be long-ranged and demand special treatment within simulations. 
A suitable tool is provided by the Ewald summation of the Rotne-Prager-Yamakawa 

mobility tensor [U [2j E3 E] . 

If the colloidal particles in a solvent are not distributed in the full 3D space, 
but rather form a thin (mono) layer, the system of interacting colloids may be 
considered quasi two-dimensional. This may be realized by either trapping the particles 
at an interface [121 EE US] or by placing them in the vicinity of a free or hard 
boundary [TJ HJ [131 Q2]> or by looking at thin fluid films only [5]. Then the question 
arises how to treat the hydrodynamic interactions in these quasi two-dimensional 
systems. In the latter case, for colloids in a thin fluid film, an experimental study 
revealed, that the two-dimensional form of the Oseen hydrodynamic tensor provides 
a suitable description of the hydrodynamic interactions of this system [5]. However, 
for colloids trapped at fluid interfaces [151 US]) or in the vicinity of interfaces as in 
the experiment described in Refs. [TJ 0] , the situation is more involved. Owing to the 
flow fields extending over half the 3D space, the system cannot be described with a 2D 
Oseen tensor with its peculiar long-ranged interactions decaying logarithmically with 
the interparticle distance. However, also the use of the full 3D Rotne-Prager-Yamakawa 
or Oseen tensors seems somewhat ill-founded, since both do not resemble solutions of 
the Stokes equation with respect to the underlying boundary conditions p3|. However, 
in the case of spherical objects these tensors still might serve as a rough approximation. 
As was advocated in ref . [T2] , the presence of a free interface separating two fluid phases, 
has only little effect on the diffusion, and thus the mobility, of spherical particles half 
immersed in both phases. The Green's function in Stokes flow (Stokeslet) for the velocity 
field for a single particle in the presence of a boundary (free interface or rigid wall) 
consists of the Oseen tensor, the free (bulk) solution, plus a mirror term [121 [331 EH] 
(method of images). Therefore, neglecting the latter term while constructing a solution 
can be considered as a leading order approximation, where the particles are assumed to 
be far from any confining boundary. This has been successfully applied to simulate the 
diffusion of particles close to an interface in experiments [TJ H]. Alternatively, the quasi 
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2D version of the mobility tensor given in Ref. [32] for particles in the vicinity of a free 
interface may be used. 

For computer simulations using Ewald summation of a 3D mobility tensor, the 
question arises how to treat the long ranged part in the third (z-) direction. Within the 
standard approach as developed by Beenakker [H], and used in the simulation studies 
presented in Refs. [H H], one has to assume periodic images also in z-direction, i.e. 
the construction of a layered system of many interfaces. Although easily implemented, 
such a procedure is computationally expensive and not necessary from a physical point 
of view. In view of the usually implemented 3D Ewald algorithm it is, however, 
unavoidable. Since we are interested in a quasi two-dimensional system which is now 
extended into the third dimension, the question arises to which extent the influence of 
the artificial periodic images in the z-direction disturbs the result of the summation. 
Indeed, one could move the periodic images in ^-direction to large distances in order to 
study their influence within the main layer of particles as function of their distance, as it 
was done for e.g. Coulomb interaction (see Ref. [19] and references therein). Although 
the effective velocities for the particles within this procedure converge to the 2D result, 
we will show in the following, that as soon as noise is considered, this approach generates 
a spurious system size dependence. Additionally, if one places the periodic images far 
away from the layer under consideration, this is even more expensive, since the sums 
within the usual Ewald formalism have to be cut off at larger K-values in reciprocal 
space. Thus for quasi 2D systems, 3D Ewald summation should be avoided. However, 
2D Ewald summation formulae have been given so far only for the Oseen tensor in an 
implicit form [20]. With regard to broader applications, note that the Oseen tensor 
suffers from not being positive definite. This renders its usage problematic as soon as 
a noise term requiring Cholesky decomposition of the tensor is present. Therefore 2D 
Ewald summation is needed for more suitable mobility tensors. Our method, as outlined 
in the following, naturally applies for the Rotne-Prager-Yamakawa tensor and also in 
the case of the quasi 2D mobility tensor of Ref. [13] . 

The paper is organized as follows. In section[2} we will first discuss the shortcomings 
of the three-dimensional Ewald summation for quasi 2D systems. Then we will formally 
derive the two-dimensional summation formula for the Rotne-Prager-Yamakawa tensor. 
As argued above, we consider this summation of the mobility tensor to be the more 
appropriate approximation than summing up the full 3D tensor since the result 
does not suffer from a spurious system-size dependence due to unphysical images 
across many additional interfaces and provides a computationally much cheaper way 
to incorporate HI within quasi two-dimensional systems. We follow the procedures 
described in Refs. [211 122], where a lower-dimensional Ewald summation has been 
developed for electrostatic and dipole interactions, and derive summation formulae for 
the aforementioned mobility tensors with periodicity assumed in two of three dimensions. 
In section [31 we demonstrate the 2D Ewald summation procedure by carrying out 
simulations of a quasi two-dimensional system. In a first step, we will show that 
conventional 3D Ewald summation leads to a divergent long time diffusion constant, as 
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the system size increases. We then apply our new summation formulae, and show, that 
the results are now independent of the system size. Finally we compare our findings 
to experimental data, and show that for the particular system under consideration, 
a reasonable agreement with the data can not be achieved using the Rotne-Prager- 
Yamakawa mobility tensor. Only upon summing the quasi 2D mobility tensor given 
by Cichocki and collaborators, simulations are found to approximate the experimental 
data. In view of the latter finding, and since the quasi 2D Ewald summation procedure 
outlined in the following could in principle be applied to any other geometrical setup or 
approximation (provided the system is quasi 2D and terms oc l/r k (k > 1) are present), 
we have derived and provide explicit formulae ready to be used in quasi-2D simulations 
for the 2D-Ewald sum of 

(i) the Rotne-Prager-Yamakawa tensor (Eqs. (I39H41I) ) 

(ii) the Oseen tensor (Appendix, Eqs. (IC.5tiC.7l) ) 

(iii) the quasi 2D mobility tensor of Cichocki et al. (Appendix, Eqs. (1C.5HC.7I) ) 

(iv) the binary Rotne-Prager tensor (Appendix, Eqs. (1D.2HD.4]) ) 

Finally, our conclusions are summarized and discussed in section HI 
2. Ewald summation for quasi two dimensional systems 

Consider a tetragonal lattice with unit cells of volume L 2 x L . The lattice is periodic in 
two dimensions (i.e. x and y) and finite in the third dimension (z). Each cell contains 
./V spherical particles, arranged to form a single layer (monolayer) parallel to the x — y 
plane. The force acting on an individual particle i will be denoted by Fj. We assume 
that no external forces are present, thus the total force on the particles in the unit cell 
cancels to zero [14] : 

N 

1=1 

If the particles are surrounded by a solvent, one expects additional hydrodynamic 
interactions between the particles. For solvents with low Reynolds number the 
motion of the colloidal particles is overdamped and inertia of the particles may be 
neglected [TT]. Concerning the implementation of hydrodynamic interactions within 
computer simulations, this leads to the use of an position-dependent mobility tensor 
for the calculation of the viscous drag of the particles [El [TT]. One particular version 
of this mobility tensor is the Rotne-Prager-Yamakawa tensor. For the integration of 
the equations of motion of the colloids, an effective velocity of each particle has to be 
calculated via 

N 

v,.:r (2) 
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with the Rotne-Prager-Yamakawa mobility tensor 



My = (67r V a)- 1 |^ 1 (1L + 

+ ^a%%l - 3fyfy) I , (^ j , r > 2a) (3a) 
M^(6^)-{(l-l^)u 



3 ry r y 

32 ary 



(i^j,r< 2a) (3b) 



M« = (67070)-% (i = j) (3c) 

and ry = |rj — rj|. The product ryfy is the outer product of the normalized vectors 
ry = Tij/rij and 11 denotes the unity matrix. Since the hydrodynamic interactions are 
long-ranged oc r^ 1 , Ewald summation has been suggested to treat the interactions of 
the periodic images [14J. This leads to a lattice sum 

N 

3=1 n 

where the second sum runs over two dimensional lattice vectors n = (n x L, n y L) with 
n/0 for ry = (indicated by the prime on the sum) and 

My(ry,n) = (677770) -1 

.'3 1 1 3 1 
X { -a- r + -a 6 -, - 11 



4 |ry + n| 2 |ry + n 

3.1 1 
" 2 a3 |r tJ + n|5 (rtJ + n)(rjJ + n) J • (5) 

Note that the definition of My for distances r < 2a (Eq. f!3bl) ). introduced to guarantee 
the positive defmiteness of My [9], is not relevant for the following. It only contributes 
to the lattice sum for n = and can be added separately. 

2.1. System size dependence of conventional 3D Ewald summation for monolayers 

Concerning the above lattice sum, we will first consider its three-dimensional analog and 
the resulting Ewald summation derived by Beenakker [H]. In order to use it, one has 
to assume periodicity in the third dimension, thus the layer of particles is reproduced 
also in ^-direction at distances n z L' [U |2]. The three-dimensional lattice sum is split 
into two sums, one in real and one in reciprocal space, respectively. We consider the 
sum in fc-space (k denotes the three-dimensional wavevector), 

1 N 

Sk = jjt £ £ M(2) (k)Fj cos(k ' r * j) (6) 

* k^O j=l 
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where the matrix M( 2 )(k) contains the /c-space part of the summed mobility tensor (see 
Eqs. (4) and (6) in ref. [H]), and concentrate on the wavevectors k with k x = k y = 0. 
Since the particles are arranged within a single layer, say at z = 0, the product k • 
vanishes for k x = k y = and all pairs of particles. Thus the cosine above equals one. If 
we move the layer in z-direction to large distances, the fundamental mode k z nin = 2ir/L z 
approaches zero. For the matrix M^(k) = M^{k z ) then follows: 



lim M (2) (A; 2 ) = lim (l-^](a- -a 3 k 3 z ) 



k\ k A z \6n ( ki 



4a 2 8a 4 / ^ * \ 4a 

t -,, ( Svra Qira 1 -> 
pa km U —r- + —j - -Qira 6 
k z ^o \ k 2 z Aa 2 3 

kk / 6ira Qua 1 o\ 

" T? (if + w ~ f a ) (7) 

Since the matrix kk/k 2 contains only one nonzero element, ((kk/k 2 ) zz = 1) the diagonal 
elements M xx and M yy of the matrix diverge as L z increases. The value of these matrix 
elements does not depend on any dynamical variable, it is constant and determined by 
the choice of system size L z and particle radius a. Therefore, after carrying out the 
3D Ewald summation, the summed mobility matrix My(rjj) consists of a dynamical 
part, and a constant part depending only on system parameters and diverging as L z 
increases (M^r^n) denotes the spatial lattice sum of the mobility matrix, c.f. Eq. 



(5) in ref. [HJ): 

M ^(^') = ^ / M 4J (r 4J ,n) 



(Gnna)- 1 (tStj + (1 - %) { E 'M (1) (r,i,n) 

^ M( 2 )(k)co S (k.r,) + ^ E M(2) W1) 



+ L 2 L 



-const. 



MJ-(rij) is in fact a 3iV x 3^ matrix, consisting of a set of 3 x 3 submatrices Dij(rij) for 
each pair of particles. After the 3D Ewald summation it contains a diverging constant 
in the xx— and yy— diagonal elements of all off-diagonal submatrices Dij,tyj{?%j)- Upon 
summing over the forces on all particles this constant part does not contribute, since 
we assumed a zero net force (see Eqs. (CQ) and fl3|). If the latter requirement is relaxed, 
i.e. in order to study sedimentation, for 3D suspensions it is necessary to include the 
backflow of the solvent as a pressure gradient in order to achieve a cancellation of 
the singular terms arising from k = [8]. However, this would not cure the system 
size dependence in this special case, since the divergence with L z does not require 
a vanishing wavevector. It is rather an artefact of the application of a summation 
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technique appropriate for three-dimensional periodic systems only. For the sum in Eq. 
([6]) the divergence with L z is therefore relevant for non-zero net forces only. However, 
the summed mobility matrix also enters the calculation of the correlated noise [TT] . 

(r i (t),r J (t + At))=2M* J (r lJ )At (9) 

thus implying a diverging width of the correlator. Within the usual Ermak algorithm, 
one has to calculate the Cholesky decomposition of M^-(r^) in order to compute 3iV 
correlated random numbers for the random displacement of the particles. Therefore, if 
M*-(i"jj) contains a diverging part, its "square root" matrix ^(r^-) with o'*jO'*T = M*j 
will also diverge and the random displacement of the particles may become arbitrarily 
large. A possible way out could be to simply subtract the divergent part, or to cut off 
the sum at a certain value, however, the latter would introduce an additional parameter 
whereas in the first case, we are not aware of a consistent argument, why the distances 
of the additional layers introduced in z— direction, should not matter. Therefore we 
consider it more appropriate to avoid this scenario by considering the system to be 
genuinely two-dimensional from the beginning. 

It is interesting to compare the setup of a monolayer of particles to another relevant 
case: particles confined between two parallel walls, rendering a three dimensional 
distribution of particles with finite extent in one dimension (see e.g. ref. [23 ] EM ] 125 ] 126] ) . 
A monolayer of colloids may be considered a limiting case of the more general 
configuration of a confined suspension of particles. However, there are also important 
differences: Within the confined geometry, particles are not only restricted to the slit by 
the walls, the different boundary conditions of the walls compared to the unbound fluid 
also alter the hydrodynamics. The distribution of particles is finite in one dimension, 
however, depending on the width of the slit, particles are able to move in the third 
direction also. The hydrodynamic interactions then also depend on the third spatial 
coordinate, and the 3D Ewald summation would therefore not diverge as in the case 
of a monolayer. The slit geometry has been dealt with in detail in ref. [23J, where the 
hydrodynamic interactions were included on the basis of a two dimensional Fourier series 
for the Green's function of the Stokes equation with the corresponding no-slip boundary 
conditions. The method has been generalized to arbitrary domains in ref. [23] as a new 
method for the computation of the hydrodynamical interactions for confined geometries, 
the so-called general geometry Ewald-like method (GGEM). The latter method relies on 
the separation of forces into local and long-ranged parts HH] [25] just as in conventional 
Ewald-summation in electrostatics. Instead of solving the Stoke equation in order to 
obtain the Green's function, our approach is rather to use an existing 3D solution 
in terms of well known bulk mobility tensors and apply them to a two-dimensional 
monolayer of particles immersed in an infinite 3D medium. In the following, we first keep 
the dependence on the non-periodic, third spatial coordinate throughout the derivation. 
Only for the final formulas, we then consider the limit of vanishing distances between the 
particles in the third direction, i.e. the spatial configuration of a monolayer. However, 
upon dropping this requirement and after some straightforward calculations, summation 
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formulae could also be obtained for a three-dimensional configuration of particles, with 
finite extent in on dimension. 

2.2. Two-dimensional Ewald summation 

In the following we will now derive the result for the Ewald summation with periodicity 
in only two of three dimensions, which will cure the above spurious dependence on the 
system size. According to the procedure described in ref. [2T], we introduce: 



n 1 J 

^ r« + n 



rij ^ 0); 



3 ' VM 3 



r« + 0); 



n 

X(r*i,£) = Yl 



exp(-i£(r i: j + n)) 



\ r ij + n l 



exp(-^(r ii + n)) 



|r»j + n| 



ftj ^ 0); 



and denote the sums for = and n ^ by $ 0) ®o an d Xo respectively: 



n 

n 

n ' ' 

x - exp(-i£(n)) 



Vo 



(«) = £ 



exp(-i£(n)) 



n 



; n ^o). 



Using these definitions, Eq. (j4j) can be written as 



N r /3 

(67rr/a)v iiCff = Fj + j ( -a$(r y ) 



a$ + 7:a 3 ^(rij) + -a 3 ^o U 



« 3 V € V a (r 4J ,^) 



*v^e (O 



4=0 



4=o 



! v^v ao (0 



(10) 

(11) 
(12) 

(13) 

(14) 

(15) 
(16) 
(17) 



4=o 
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where the denotes the gradient with respect to £ and is evaluated for £ = 
Note that the gradients appear with the opposite sign to compensate for the i 2 due to 
differentiation of Eq. (TT2D and Eq. 1115]) . 

The sums appearing in Eqs. ( |T0|) -( TT3|) will now each be transformed into two rapid 
converging sums in real and reciprocal space. We repeat the formalism described in 
Ref. [21], in detail for the first component $(iy,) and give the results for the other parts. 
For this purpose we use the definition of the Gamma function 



1 1 

^ " T>) Jo 



l dt (19) 



and Poisson's summation formula in two dimensions for a sum of Gaussians 

71 K 2 

^exp(-|p + n| 2 t) = _^ eX p(zKp)exp(- — ) (20) 

n K 

where p = (x,y) and K = (k x ,k y ) are two-dimensional vectors in real and reciprocal 
space respectively. Note that the three-dimensional version of this summation formula 
differs only be an additional factor of y/rr/Lt' 1 ^ 2 on the right hand side, stemming 
from the underlying Fourier transformation of the sum of Gaussians. Using the 3D 
summation formula instead, e.g. in order to derive the original result of Beenakker for 
the 3D Rotne-Prager tensor [T3], would therefore not change the general outline of the 
calculation, but lead to different integrals in the following. 
Inserting Eq. (TIH} for s = 1/2 into Eq. (TH)]) leads to [21] 



oo 



$ = E f(l72) I rl/2ex PH^' + n\ 2 t)dt (21) 

The integral may be split up, by introducing a convergence factor a [21] (still to be 
determined) 

1 f°° 

$ = ) —= / t~ 1/2 exp(-|rij + n\ 2 t)dt 
„ V 71 " JoP- 

i r a2 

+ V -= / t~ 1/2 exp(-|i-y + n\ 2 t)dt (22) 
The first integral can be evaluated, for the second we apply the summation formula Eq. 



n 



erfc(a 


ry + n|) 




r i.i 


+ n 





+ T t~ 3/2 J2 ex P (iK Pii ) exp(-^ - z 2 3 t)dt (23) 

K 

where erfc(x) = 1 — erf(x) is the complementary error function and = (pij, Zy). The 
integral will contain a singularity for K = 0. This term deserves a special treatment, 
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therfore, we separate the term with A = from the sum to isolate the singularity 
As will be seen later on we can use the absence of external forces (Eq. (CQ)) to get rid 
off all singular and constant terms [HJ I2T] . 



\ r ij + n l 



L 2 



The first integral evaluates to 



(24) 



i>a 2 

/ t~ 3/2 exp(-zlt)dt 
Jo 



2v / 7r-2iierf (aza) e a Z{ i + lim 

a t^o+ y/t 



(25) 



while the second integral can be performed using the substitution u 2 = 1/t [21]: 

pa 2 

J t-^ 2 exp(-—-z 2 J t)dt = 



A 



-Kzi 



A 



erfc(^ — azij) + e Zij erfc(^- + az 



K 



2a 



2a 



(26) 



Putting all the pieces together, one ends up with 

erfc(a|rjj + n|) 



_7T_ x - exp(zKpij) 
i 2 ^ 



A' 



A 



+ e KziJ erfc(— + az {j ) 



2y^F 

L 2 



V^z^erf (azj 



1 _ Q 2 2 2 e ~« 

— e « — hm 
a 



(27) 



Accordingly, for $o (Eq. ( 1141) ) one finds 



erfc(a 


|n|) 




n 





A 



20F + 20F 1: 



E r rfc( 2^ 

2a 



L 2 a L 2 t^o+ y/t y/n 



(28) 



where the term for n = had to be inserted in the sum of the integral on [0, a 2 } and 
subtracted separately in order to apply the Poisson summation formula [2T] . 
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For the terms ^(r^) and \l/o the procedure is similar, so we just give the results. 
Note, that no singular terms will appear while performing the required integrations. 
Use Eq. flEJ} with s = 3/2 leads to 

*(r y ) = — V / t 1 / 2 exp(-|r ii + n\ 2 t)dt 
V 71 „ Jo? 

+ ^E ex P( zK ^)/ r 1/2 exp(-^-zJt)rft (29) 



and 



2 Z" 00 
= — 5 / t 1/2 exp(-|n| 2 t)dt 

V 71 " „ Ja 2 



K ^° 



^ 2 Vio " V 4t'- 3v^F 



(30) 



The G and x terms require a different form of the Poisson summation formula |21j : 
^exp( - |p + n| 2 t-i£(p + n)) 

n 

= £ ex P( iK P) ex P(- ^ K ^"^ ) ( 31 ) 

K 

which takes the additional ^-dependence into account. Using Eq. (TT9]) with s = 3/2, 
Eq. (I12p can be written in the following form: 

2 x - 

= -f= V exp(-?^(r ij - + n)) 
V 71 " 



n 



X 



/•oo 

/ t 1/2 exp(-|i\,- + n\H)dt (32) 



Again, one splits up the integral and applies Eq. (131)) in the second integral. Thus 

2 x 



n 



x / t 1/2 exp(-|r ij + n| 2 t)dt 

K 



^ t- 1/2 exp(- |JV ^' -4-t)dt (33) 
where £ p denotes the two-dimensional x- and y-part of and accordingly for 0o(O 

c\ POO 

e (0 = ^Vexp(-^(n)) / t 1 /2 exp( _| n |2 t) ^ 

V 7 ^ „ Ja 2 

r 2 i/ 2 / + 4" 3 / 



|K + ^| 2 

4t 



L 2 ^Jo ^ 4t 3^ 
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The calculation for xi r iji£) is carried out with s = 5/2 and after some short 
manipulations one finds 



4 x - 

x( r y > £) = g-y= 2^ ex PH£( r t> + n )) 

* n 

poo 

x / t 3/2 exp(-|r ij + n\H)dt 

Jo? 

K 



X 



12 



t 1/2 exp(- ' -4-t)dt (35) 



4 T 00 
Xo(0 = Vexp(-z£(n)) / t 3/2 exp(-\n\H)dt 

^ 1/2 , l K + CI\, 8« 3 , x 



K 

The calculation of the components for = 0, i.e. $o, ^o, @o and Xo (see Eqs. 
(I14til7p ) is only needed to extract the constant term arising from the constraint K ^ 0. 
This term only contributes for % = j, since it corresponds to the additional part summed 
up for = and n = 0. Therefore it is placed outside the sum over j. The remaining 
part can be absorbed into the main formulae by allowing = 0. Note that for and 
X this constant term does not contribute, since it doesn't depend on £ and only the 
matrix elements with gradients are used. A similar singularity as in the integration (Eq. 
( 125]) ) will appear when integrating V^V^O (see also Eq. ( 1A.1I) in the appendix). 

We are now left with expressions for $, G and \ as sums in real and reciprocal 
space. It remains to evaluate the gradients with respect to perform the remaining 
integrals and collect all parts from Eqs. (I27H30P and Eqs. (I33H36P in order to insert 
them into Eq. (j!8j) . We list all components and integrals in the appendix. To obtain 
the final formula, we take the limit z^ — > for all pairs of particles, i.e. all particles 
remain in a plane parallel to the (x — y)-plane. We introduce the following definitions: 

i^ n =|r*; + n|, R jn = - , K = § . (37) 

\r.ij + n| K 

Since we are only interested in the limit = 0, we can assume without loss of generality 
that z = for all particles, i.e. = pij, and thus consider Rj n , Vj C ff and F, to be two- 
dimensional only from now on. Before writing down the sum of all components we note 
that several terms cancel and, additionally, all divergent and constant parts stemming 
from Eqs. (I27H30I) and (I33H36I) do not contribute to the final result. In the case of a net 
zero force on the particles in the system, (see Eq. (CQ)), any product with a (diverging) 
constant also vanishes. As discussed in Ref. [8], the cancellation of the singular terms 
stemming form K — can be achieved, even in the case when the average force on the 
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particles is not zero. To this end, one has to introduce a pressure gradient representing 
the backflow of the fluid. The relevant physical quantities are then the velocities relative 
to the backflow of the fluid |8J. Taking all similar parts together, one then arrives at a 
rather compact result: 

aa f 3 2 2 2 \ 
(67r?7a)v iiCff = Fj - —j= I - + -a a I F, 



N ( 

+ E E 

j=l I n 



. / 3 erfc(aiL n ) 
' -a — 

4 R jn 



+ 



, _ fv 2p2 



1 3 erfc(a-Rj n ) 



11 + R,j n R,j 

11 3I^/j jilv.jjj 



aa 



7i 



2a 2 a 2 



+ ^cos(Kr,,) (l^erfc^J 

30-^/7? 



]L-(i-ia 2 iT 2 )KK 
v 2 3 ; 



2L 2 a 

If we cast this result into a similar form as in Ref. 

A' 



J 



KKJ | F 

the final formula reads 



(3* 



(67070) Vi , eS = ,Af(1) (%n)Pi 



j=l n 



^^M( 2 )(K)cos(Kr 4 ,)F, 
i=i k^o 



(39) 



The prime in the first sum indicates, that for r,y = the terms with n = are omitted. 
We used the definitions: 



M«(i 



11 



3 _i 1 o o 

-ar H — a r 

4 2 



J erfc(ar) 



3 _2 -1/2 / 2 2\") 

a w 7r 1 exp(— a r )) 



rr 



3 1 3 O Q 

-ar a r 



J erfc(ar) + (— 3a 3 ar 2 



-aa — 2a 3 a 3 ) it i/z exp(— a z r 



-V2, 



(40) 
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for the real part of the tensor and 



M (2) (k) = fl 1 2a erfc 

«~ (} 2 

— kk < [a a 

H 3 

+ aa _1 /c7r _1 ' 2 e: 




(41) 



for the summation in Fourier space. This completes the derivation of the quasi two- 
dimensional Ewald sum of the Rotne-Prager-Yamakawa mobility tensor. 

The main advantage of the Rotne-Prager-Yamakawa tensor compared to the Oseen 
tensor is its positive definiteness [9], needed for the Cholesky decomposition as used in 
Brownian dynamics simulations [IT]. Therefore, one also has to consider the part of the 
Rotne-Prager-Yamakawa tensor for r^ < 2a (Eq. (23) in Ref. [9]). However, this term 
does not contribute to the long-ranged part of the Hydrodynamic Interaction and thus 
can be simply taken out from the real space lattice sum for n = in Eq. (j38p or Eq. 
(139]) and added separately. For the sum over the periodic images, this term will not 
appear since jr^- + n| > 2a is always fulfilled. 

3. Results from simulations 

In order to test the Ewald summation in a quasi-2D system using the procedure outlined 
above, we performed Brownian dynamics simulations with a single layer of N colloids 
(Radius a) within a fluid phase in a box with side lengths L x = L y = L, L z . This 
system resembles a model for the setup used in the experiments of refs. [TJ [17] where 
paramagnetic colloidal particles were placed atop a flat and stabilized air-water interface 
of a suspended droplet. The colloids are fully immersed in the fluid phase and, due to 
gravity pulling them downwards, just stay in the vicinity of the interface [27] without 
perturbing the latter considerably. Thus the colloidal particles constitute a quasi two- 
dimensional system. More details can be found in refs. [HUT]. In the model, the layer of 
particles is placed parallel to the (x — y)-plane at z = L z /2, with only in-plane motion 
allowed. The colloids interact through a repulsive potential Vd/k^T = T/d 3 , where d 
is the distance between each two particles scaled by the mean interparticle separation: 
d — r l\/~Q {q denotes the 2D number density of the colloids). In the experimental setup, 
this dipole repulsion is generated and controlled by an external magnetic field applied 
perpendicular to the layer of colloids. 

First, we validate our simulation by attempting to reproduce the Brownian 
dynamics data published by Rinn et al. pp. We therefore set Y = 8.2 and a = 2.35/im for 
a system of N = 100 colloids at a number density of g = 3.24 • 10~ 3 , which corresponds 
to an area fraction of rj = 0.056. The time is measured in units of r = 1/(qDq), where 
D is the single particle short-time diffusion constant, extracted from the simulations 
by extrapolating the mean-squared displacement towards t — > [I] . In our simulations 
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we used the bulk value for the diffusion constant (D = (6^770) _1 ) . If the scaling of 
the mean-squared displacement with D holds, the extrapolated short-time diffusion 
constant Dq should agree with this value. Indeed, for Brownian dynamics simulations 
we find Dq = D and the data agrees with ref. [1] (see Fig. [1]). Furthermore, we 
compare in Fig. [1] the mean-squared displacement (scaled by 4tDo) for simulations 
with varying system size L z . First it should be noted that the the extracted short 
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Figure 1. Comparison of the scaled mean-squared displacement for Brownian 
dynamics simulations of N = 100 colloidal particles (radius a — 2.35/iin, 2D 
number density g = 3.24 • 10 -3 /zm -2 ) with a repulsion strength of V = 8.2 without 
hydrodynamical interactions (open circles) and with HI included via 3D Ewald 
summation for a varying longitudinal system size L z . Experimental data and additional 
BD data (line) were taken from ref. [T] . Errorbars are of the order of the corresponding 
symbol size and have been omitted for clarity. A simulation with L z w 7L (not 
shown) would accidently match the data from the experiment. Such a fortuitous choice 
could well be the underlying reason for the good agreement between simulations with 
hydrodynamic interactions and the experimental data, as reported in ref. [T]. 



time diffusion constant D differs from the corresponding bulk value D, if 3D Ewald 
summation is applied. The extracted values increase with increasing system size in 
z— direction. Additionally, as can be seen from the simulation data presented in Fig. [TJ 
that the mean squared displacement, scaled by the corresponding short time diffusion 
constant Do, also increases with increasing system size. Thus the simulations confirm 
the finding of a diverging mobility matrix for the 3D Ewald summation of this system. 
Neglecting noise, the 3D summation method remains valid. If we then compare the 
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effective velocities of single particles, we found that for our current setup, the distance 
of the periodic images in z-direction had to be scaled by a factor of two relative to the 
original size L z = L X;V in order to reduce the impact on the single layer and to converge 
to the result from 2D Ewald summation according to Eq. (|38l) . Note, that for large 
systems, such a stretching of the third dimension could be unnecessary, since the size of 
a cubed box could be already sufficient. However, this has to be checked for each setup 
individually. 

As a second test, we perform simulations with hydrodynamical interactions 
included, but without applying any Ewald summation procedure. We extract the long- 
time self diffusion constant Dl by running the simulation for much longer times t = 1.2r 
and fitting the scaled mean-squared displacement for times t > 0.95r to a constant. Fig. 
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Figure 2. Long-time self diffusion constant (scaled by Do) extracted from simulations 
with hydrodynamical interactions, with and without 2D Ewald summation. As the 
number of particles in the system increases, Dl approaches the value as obtained from 
simulations with 2D Ewald summation. Lines arc drawn to guide the eyes, the dashed 
line corresponds to a polynomial fit. Error bars correspond to the statistical error 
obtained from averaging over many simulation runs. 



|2] depicts the long-time diffusion constant (scaled by Do) for simulations with the same 
parameter setup as before and for various system sizes and constant number density. As 
the number of particles increases, for the Brownian dynamic simulation with HI, Dl is 
found to decrease and approach the limit set by using 2D Ewald summation. The system 
size dependence of D L without Ewald summation becomes clearly visible. In order to 
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roughly characterize the convergence, we extrapolated this decrease using a polynomial 
fit for < N < 700. It turned out that a system containing N > 25000 particles would 
yield a similar result for Dl/Dq as obtained in simulations with 2D Ewald summation. 
For the latter ones, we do not observe any significant dependence of on the system 
size. 

Of course it is now tempting to compare the 2D Ewald summation results to the 
experimental data from ref. pp. However, one has to keep in mind, that the situation in 
the experiment is different. The particles are adsorbed to a free interface, whereas there 
is none in the simulations. The authors suggested to increase the hydrodynamical radius 
in the simulations, that is, the value of the particle radius within the calculation of the 
mobility tensor was increased by a factor of two. Fig. [3] depicts the scaled mean-squared 
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Figure 3. Comparison of the simulated scaled mean-squared displacement for the 
same setup as in fig. [T] to experimental data as obtained by the Experiment of Rinn 
et al. jTj. Only representative error bars are shown for simulations. The error bars 
for the Brownian dynamics simulation without hydrodynamical interactions (stars), 
are smaller than the symbol size. The simulation data stemming from BD simulations 
with hydrodynamical interactions and 2D Ewald summation of the quasi 2D mobility 
tensor of Cichocki et al. [13] (triangles) has been scaled by their different short time 
diffusion constant Dq = 1.38Z?. 

displacement for the 2D Ewald summation method and the experimental data of ref. pQ . 
The data cannot be reproduced with this actual implementation. We refrained from 
trying out various hydrodynamical radii, since this would introduce a free parameter 
to the system. If we use the mobility tensor by Cichocki et al. [T5] , derived from the 
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two-sphere mobility tensor for particles close to a free interface as an asymptotic series 
in powers of 1/R (c. f. Appendix C Eqs. fIC.lal) and fIC.lbj) ). the simulation data 
overestimate the experimental values. However, the short time diffusion constant -Do 
extracted from the simulation is also increased, compared to the simulations based on 
the Rotne-Prager mobility tensor. If we normalize the quasi 2D simulation data with 
respect to the extracted value for Do, we find a rather good agreement. For longer 
times, however, the simulation data starts to deviate slightly. Note that the deviation 
of the extrapolated short time diffusion constant D compared to D = (Gnrja)" 1 may be 
anticipated by inspection of Eqs. (IC.lal) and (IC.lbj) in the appendix, since the matrix 
Qi for the self-diffusion already deviates by a factor of ~ 1.38 from the unity matrix in 
the Rotne-Prager case. Therefore, this deviation is only present in the quasi 2D case, 
for all other simulations the value of Dq agrees with the diffusion constant D initially 
plugged in. 

Concerning the computational cost of the simulations, we found that using the 
2D Ewald summation method lead to a reduction of at least an order of magnitude, 
compared to its 3D version, and depending on the particular choice of the convergence 
parameter a. For a typical choice of a = 2/L the gain in speed was around a factor 
~ 15. Note that the computational cost of the Ewald summation may vary, depending 
on the choice of a. 



4. Summary and conclusions 

In summary, we have provided Ewald summation formulae for quasi two-dimensional 
systems for the Rotne-Prager- Yamakawa mobility tensor and variants. We 
demonstrated, that for quasi-two-dimensional systems, 3D Ewald summation leads to 
a spurious system size dependence, stemming from the summation in the direction 
perpendicular to the 2D layer of particles. This problem was solved by summing 
in two dimensions only, and, additionally, using the resulting formulae to calculate 
hydrodynamic interactions in computer simulations of quasi 2D systems was found to 
be much more efficient, due to the avoidance of summation in the third direction. 
We further found that the asymptotic value of the long time diffusion constant for 
large systems could already be obtained for rather small systems using the 2D Ewald 
summation procedure. We demonstrated, that the 2D Ewald sum of the quasi- 
two-dimensional analog of the Rotne-Prager mobility tensor given by Cichocki and 
collaborators may be used to reproduce experimental data quite well . Together 
with recent advances for an approximate and efficient treatment of HI in computer 
simulations [28], inclusion of HI and the proper treatment of their long-ranged 
characteristic becomes feasible even for large systems in quasi 2D simulations of colloidal 
suspensions. 

J.B. thanks M. Oettel for fruitful discussions and the German Research Foundation 
(DFG) for the financial support through the Collaborative Research Center (SFB-TR6) 
"Colloids in External Fields" Project N01. 
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Appendix A. Gradient terms and Integrals 



Here we list the evaluation of the gradients V^V^ for the integrals appearing in the 
Equations for and % (Eqs. (I33H36I) ). as well as the resulting integrals that need to be 
performed in order to obtain Eq. (138]) : 
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The integral containing t 3 / 2 has already been calculated (see Eq. (126ft ) . the remaining 
integrals can also be performed: 
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Taking the limit z — > is straightforward for all terms except for the integral in Eq. 
(1A.3j) . This evaluates to: 
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Appendix B. Ewald sum of the Oseen tensor 

For the Ewald sum according to Eq. (jlj) of the Oseen tensor 



O iJ = (8 7 rr / )- 1 r^(ll + f^), (i ^ j) 
On = {Qirria)~ l \ } (i = j) 



(B.la) 
(B.lb) 



instead of Mjj, we neglect all terms in Eq. (!38|) involving a 3 , since these are the 
additional terms of the Rotne-Prager-Yamakawa tensor compared to the Oseen tensor. 
This leads to: 
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Which also may be casted in a similar form as Eqs. ( l39ll4Tj) : 
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for the real part of the tensor and 



M( 2) (k) = u{2aerfc(A^| 



3tt 



— kk s 

for the Fourier space. 
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Appendix C. Ewald sum for the quasi 2D mobility tensor of Cichocki et al. 



The mobility tensor reported in Ref. [13] for a quasi 2D system of spherical particles 
close to a fluid interface differs from the 3D Rotne-Prager-Yamakawa tensor by a factor 
of two for the terms proportional to l/r^ and by a factor of q = 5.59027 for terms 
oc l/rfj. Additionally, it uses different matrices for the self mobility and for terms 
oc 1 / rfj instead of the unity matrix: 



M, l3 = (6n V a)- 1 ||or^(U- 



r ij r ij; 



M u = (6^)"^!, (i = ]) 



(i * h 



with q = 5.59027 and matrices Qi and Q 3 [13] 
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-0.319658 



(C.3) 



reflecting the presence of a free interface (the additional scaling factor q and the matrix 
Q3 may be derived from Eq.(8) in Ref. [13] by casting it into the Rotne-Prager form 
(Eq. fl3al) )). As a consequence of that, the resulting summation formula involves more 
terms, since the different matrices avert cancellations. The corresponding analog to 
Eq. fl38l then reads: 
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and accordingly 
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with the definitions: 

M%(r) = 1 [ (l^ 1 + ^ 3 ^ 3 Q 3 ) erfc(ar) 
+ ga 3 ar" 2 7r" 1/2 exp(-a 2 r 2 )Q 3 } 
+ f f | r^ ar_1 — ^? a3r ~ 3 ^ erfc(ar) + (— 3qa 3 ar~ 2 
+ 3aa - 2gaV) 7r~ 1/2 exp(-a 2 r 2 ) } (C.6) 
for the real part of the tensor and 

<i>M = i { (20 + i 9 A- 2 (U - Q 3 )) erfc 

-§^W'/>exp(-il) (U-Q 3 )}^ 

+ "^- 1/2 ^(-3}ll < C7 > 

in Fourier space. 

Appendix D. Ewald sum for the binary Rotne-Prager Tensor 

For binary mixtures of particles with radii aj G {do, ai}, we replace the particles radius 
a by a,, and within the sum over particles with radius a J; each factor a 3 is replaced by 
f(a 2 + a 2 ) [21: 

My = (evrr/ai)" 1 ||Oi'"y 1 (lL + fyfy) 

+ |(a 2 + r,j)r,/n - 3fyfy)} , (i ^ j) (D.la) 
Mii= (evrr/a,)- 1 !!, (i = 3) (D.lb) 
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Carrying out this procedure for the previous results (Eqs. (I39H41P leads to: 
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With the corresponding definitions: 
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for the real part of the tensor and 
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for the summation in Fourier space. 
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